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In recent years, ultracold atoms in optical lattices have proven their great value as quantum simulators for 
studying strongly correlated phases and complex phenomena in solid-state systems. Here we reveal their po¬ 
tential as quantum simulators for molecular physics and propose a technique to image the three-dimensional 
molecular orbitals with high resolution. The outstanding tunability of ultracold atoms in terms of potential 
and interaction offer fully adjustable model systems for gaining deep insight into the electronic structure of 
molecules. We study the orbitals of an artificial benzene molecule and discuss the effect of tunable interactions 
in its conjugated n electron system with special regard to localization and spin order. The dynamical time scales 
of ultracold atom simulators are on the order of milliseconds, which allows for the time-resolved monitoring 
of a broad range of dynamical processes. As an example, we compute the hole dynamics in the conjugated tt 
system of the artificial benzene molecule. 


The structure of molecules is usually determined by x-ray 
or electron diffraction. Current advances with femtosecond 
pulses allow for the time-resolved observation of the atomic 
positions [1]. In general, orbital wave functions are much 
harder to image. For simple molecules, the reconstruction of 
the HOMO has been achieved recently by electron momen¬ 
tum spectroscopy [2], by using laser-induced electron diffrac¬ 
tion [3], and by higher-order harmonic generation [4, 5]. For 
hydrogen atoms the nodal lines of Stark states were resolved 
via photoionization microscopy [6]. State-of-the-art scanning 
probe microscopy allows resolving the HOMO as well as 
identifying the chemical bonds [7] for benzene-based com¬ 
pounds attached to surfaces. 

Experimental access to the electronic structure of 
molecules and their dynamics is essential because even for 
relatively small molecules the full many-particle problem is 
not computable using classical computers. This has brought 
forward the idea of quantum computation for molecules [8, 9]. 
The quantum computation of a hydrogen molecule in a min¬ 
imal basis has already succeeded [10, 11], but estimates for 
more complex molecules are not promising. While the num¬ 
ber of required qubits appears manageable, the estimated 
number of quantum gate operations (10 18 for Fe 2 S 2 ) is many 
orders of magnitude larger than currently available [12]. 

In the past decade, ultracold atoms in optical lattices have 
been established as quantum simulators for condensed mat¬ 
ter systems [13, 14]. In addition, ultracold atoms were pro¬ 
posed as simulators for different systems such as neutron stars 
[15], black holes [16, 17], quarks [18], or atoms [19]. Re¬ 
cent experimental developments include prospects such as the 
single-lattice-site imaging and addressing [20, 21] as well as 
the deterministic preparation and detection of few-atom sam¬ 
ples [22, 23]. 

Here we show that ultracold atoms can be employed as a 
quantum simulator for molecules using existing experimen¬ 
tal techniques. In this setting, the ultracold atoms mimic the 


electrons in a molecule, whereas the optical trapping poten¬ 
tial takes the role of the nuclei. We demonstrate that ultra¬ 
cold atoms can serve as a tunable model system allowing the 
investigation of open questions in molecular physics. In con¬ 
trast to the quantum computation approach aiming at the exact 
calculation of molecular energies, the present quantum emu¬ 
lation approach uses model systems to investigate specific in¬ 
teraction effects and dynamical processes in moleculelike sys¬ 
tems. As a concrete example, we focus on a model system for 
benzene and compute its molecular orbitals for vanishing in¬ 
teraction. For the nonsolvable interacting problem, ultracold 
atoms can serve as a quantum simulator for static and dynam¬ 
ical electronic properties in molecules. We demonstrate nu¬ 
merically that the conjugated i r electron system shows strong- 
correlation effects such as localization and spin order, even 
when neglecting the interaction with inner particles. On the 
basis of this subsystem, we also reveal that ultracold atom 
simulators promise unique insight into electronic femtosecond 
dynamics. We show that momentum mapping in combination 
with phase retrieval allows the imaging of molecular orbitals 
with 1-2 orders of magnitude better than the intramolecular 
distances. We explain how to use the outstanding tunability of 
interaction and potential for studying electronic interactions 
and the dynamics in artificial molecules. 

I. CREATING ARTIFICIAL MOLECULES 

As an example, we discuss how to create an artificial ben¬ 
zene molecule and compute the orbital wave functions. In real 
benzene, the molecular structure is formed by a ring of six 
carbon atoms and six hydrogen atoms (Fig. lc). Thereby, the 
molecular symmetry is essential for the formation of molecu¬ 
lar orbitals, where benzene belongs to the point group D^. In 
our case, the idea is to simulate the electrons in a molecular 
structure via ultracold atoms in a tailored optical dipole poten- 
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FIG. 1. Molecular orbitals of artificial benzene, (a) Illustration of a honeycomb optical lattice superposed with a dipole trap. The hexagonal 
ring and the adjacent sites form a trapping potential for an artificial, benzenelike molecule, (b) Calculated molecular single-particle orbitals 
of the artificial benzene molecule with low-lying s orbitals (1-6) on the inner ”C” ring, hybridized sp 2 orbitals (7-18,25-30) including the 
adjacent ”H” atoms, and p z orbitals (19-24) forming a conjugated tv system. The orbitals are plotted for the lattice depths Vo = IITTr and 
Vo* = 35 Er in units of the recoil energy Er of the lattice wavelength. The isosurfaces depict the orbital wave functions at 0.4/av^ (orange) 
and —0.4/a(green) with the lattice constant a of the honeycomb and b of the orthogonal lattice, (c) Chemical structure of the benzene 
molecule CeH6. 


tial that mimics the Coulomb interaction with the nuclei. The 
optical potential is imposed by the ac Stark shift of a laser field 
and can be created for the example of benzene by superposing 
a hexagonal optical lattice VoVh c (±, y ) [24] with a tight dipole 
trap Vdi p (x, y) as shown in Fig. la (see Appendix). A similar 
potential can, e.g., also be created by adding a triangular su¬ 
perlattice to the honeycomb lattice (a ’’benzene lattice” [25]) 
or by using a spatial light modulator, which renders the possi¬ 
bility of almost arbitrary two-dimensional potentials [26-28]. 
In the orthogonal direction a one-dimensional lattice Vm with 
depth Voz or a light sheet is applied. A controlled number of 
fermionic atoms (e.g. 6 Li atoms with spin states ±1/2) can 
be loaded into this optical potential [22] (see Appendix). 

The molecular orbitals of the artificial molecule are com¬ 
puted using the plane-wave expansion method for vanishing 
interaction (see Appendix for technical details and param¬ 
eters). The resulting single-particle orbitals are shown in 
Fig. lb. The six lowest s orbitals represent an energetically 
separated shell within the inner ”C” ring, whereas the higher 
orbitals (7-18, 25-30) are sp 2 -hybridized. The orbitals (7-11) 
and the antibonding orbital (13) form p orbitals pointing along 
the C ring. The outward-pointing p orbitals (12) and (14-18) 
lower their energy by maximizing the overlap to the s orbitals 
of the H atoms. The hybridized orbitals (25-30) show a nodal 
plane between H and C ring and thus have a higher energy. 
Energetically in between lie the delocalized p z orbitals (19- 
24) with a nodal plane at z = 0. Benzene is one of the prime 


examples for this type of delocalized conjugation of a tv elec¬ 
tron system stabilizing the planar ring structure. This property 
is commonly referred to as aromaticity. Accounting for in to¬ 
tal 42 electrons, in benzene three of the six p z orbitals are 
occupied. We can tune the energetic order of the orbitals, by 
tuning the parameters of the potential, e.g. strength and width 
of the dipole trap as well as the lattice depths Vo and Vo z (see 
Fig. 3a). 


II. IMAGING OF MOLECULAR WAVE FUNCTIONS 

Quantum gas microscopes have demonstrated optical imag¬ 
ing with single-site resolution in an optical lattice [20, 21]. 
However, the imaging of molecular wave functions of artifi¬ 
cially created molecules requires an optical resolution several 
times better than the optical lattice spacing. The idea is to 
overcome the limited spatial resolution by imaging the par¬ 
ticles in momentum space. In ultracold atom experiments, 
this can be achieved by turning off the optical potential, let¬ 
ting the particles freely expand and imaging the particle den¬ 
sity /?tof(i*) after a certain time of flight. In the language of 
molecular physics, this would be equivalent to a sudden re¬ 
moval of all nuclei. During the expansion, the interactions 
quickly become negligible, because there are only few par¬ 
ticles with high momenta. This allows for a free expansion 
in two dimensions with an unaltered lattice in the z direction 
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avoiding problems with a limited depth of focus of the imag¬ 
ing system. In addition, the interaction among the particles 
can be tuned to zero via a Feshbach resonance during the ex¬ 
pansion. For suitably long expansion times t (see Supplemen¬ 
tal Material [25]), the momentum density can be expressed as 
p p ( q) oc Ptof (pf/ tyi) with momenta p = ftq/a, a the lattice 
constant of the honeycomb lattice, and h the Planck constant. 

In the example shown in Fig. 2b, we assume that we can 
image the momentum density after the time of flight with 
P 2 = 49 x 49 discrete momenta using a so-called pinning 
lattice for imaging (see Appendix). Initially preparing one 
particle (or two with opposite spin without interaction) only 
the lowest orbital n = 1 is occupied and per experimental 
sequence only 1-2 momenta are recorded. To show the fea¬ 
sibility of the proposed imaging method, we use a random 
number generator with an accumulated number of 1000 mo¬ 
menta. Note that in principle satisfying results with fewer ob¬ 
served momenta are possible (see Supplemental Material [25], 
Fig. S4). Figure 2b (first row) shows that the statics is clearly 
sufficient to map out the momentum distribution of the or¬ 
bitals. The phase of the distribution can be retrieved by iden¬ 
tifying the nodal lines of the momentum density (p p (q) = 0) 
and mapping it piecewise to the momentum wave function 
with p(q) = ±^/p p (q) (second row). Since the potential 
is symmetric, p( q) can be assumed to be real. The compo¬ 
nents of p(q) are the Fourier coefficients of the wave function 
and therefore allow the direct transformation to orbital wave 
functions in real space via 


ip(x,y) = X]p(q) |q) (l) 

q 


with plane waves 
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i 
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lattice constant a and the number of grid points S per recip¬ 
rocal lattice vector. The resulting density is plotted in the 
third row of Fig. 2b. The figure demonstrates that the effec¬ 
tive resolution of the imaging method is very high. In prin¬ 
ciple, the resolution is given by the diameter in momentum 
space. However, since the momentum wave function drops 
exponentially to zero for higher values of q, we can assume 
p( q) ~ 0 for |q| > 3.5. Since the diameter in momentum 
space is not limiting, the spatial resolution in the experiment 
is mainly influenced by the statistics of the momentum den¬ 
sity (see Supplemental Material [25], Fig. S4). In contrary, 
the maximal diameter of the artificial molecule to be imaged 
is given by the resolution in momentum space and estimated 
as P/ (2g™ ax — 2 q™ m ) using the Nyquist theorem (3.5a for our 
example). 

In analogy, higher single-particle orbital wave functions can 
be imaged as depicted in Fig. 2c. The preparation of higher 
orbitals is analogous to the preparation of higher bands in op¬ 
tical lattices via Bragg transfer [29, 30] or via amplitude mod- 
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FIG. 2. Imaging of molecular wave functions via momentum-space 
mapping. The first row shows the momentum density in two dimen¬ 
sions, the second row the momentum wave function after phase re¬ 
trieval, the third row the density (the sign of the wave function is in¬ 
dicated), and the fourth row the isosurface of the three-dimensional 
wave function for ±0A5/aVb. (a) Theoretical calculation for the 
lowest orbital (n = 1) of the artificial benzene molecule, (b) Simu¬ 
lated time-of-flight measurements with 1000 observed momenta in 
the x-y plane (randomly distributed) for the orbitals n— 1 and (c) 
n— 12. (d) The momentum density in y-z direction with 100 random 
momenta for the lowest p z orbital (n = 19), integrated momentum 
wave function and density in 2 direction. In combination with the 
reconstructed wave function in x-y coordinates (see (b)), this allows 
one to retrieve the three-dimensional wave function for the lowest p z 
orbital (analogous for (b) and (c) using the lowest Wannier orbital in 
z direction). 


ulation [31]. Alternatively, increasing the number of parti¬ 
cles in the artificial molecule automatically occupies higher 
orbitals. In this case, the recorded images correspond to the 
sum of the respective momentum densities (for vanishing in¬ 
teractions). Imaging the wave functions in three dimensions 
requires measuring the momentum distribution also in the z 
direction, where Fig. 2d depicts the reconstruction of the p z 
orbital. Here, a statistics of 100 particles is sufficient to re¬ 
solve the orbital structure. In our case, the separable den¬ 
sities in the x-y and z direction can be simply combined to 
the three-dimensional orbital wave functions shown in Fig. 2 
(fourth row). 

























4 


III. QUANTUM SIMULATION OF 
THE MANY-BODY PROBLEM 


The single-particle energies of the orbitals are plotted in 
Figs. 3a and b, where the six p z orbitals are indicated in red. 
Because of the separability of the potential in the z direction, 
these orbitals are analogous to the six lowest s orbitals but 
with a nodal plane for 2 = 0. The lowest and highest orbital 
are symmetric (A 2u ) and antisymmetric (BJ g ) under rotation 
(compare Fig. lb), whereas the second-third and fourth-fifth 
orbitals are doubly degenerate and labeled as Ei g and E^. The 
lower three orbitals have a bonding character and the higher 
three orbitals an antibonding character. While molecular po¬ 
tential curves are usually plotted against the nuclei distance 
Ro , the distances are fixed in an optical lattice potential. For 
the same reason the repulsion of the nuclei, dominating the 
potential curves at Ro —>> 0, is not present in optical potentials. 
However, even though the positions are fixed in our case, we 
can effectively change the intramolecular distances by tuning 
the depth Vo of the lattice potential. A deeper lattice potential 
decreases the density between the sites and is analogous to a 
larger spacing R 0 between the nuclei coordinates. 

A quantum simulator for molecular problems is needed for 
non-vanishing interactions, since solving the many-particle 
problem using the full configuration-interaction method is 
not feasible for complex molecules even accounting for the 
molecular symmetries. Unlike electrons, neutral atoms do not 
interact via Coulomb force, but by either a short-range contact 
potential or a long-range dipolar potential [32]. For ultracold 
atoms, the interaction strength between ultracold atoms can 
be tuned in a wide range via an external magnetic field due to 
the existence of Feshbach resonances. As a concrete example, 
the 8-wave scattering length a s for 6 Li atoms can be tuned be¬ 
tween —0.26 and 0.46, where 6 = A/2 = 532 nm is the lattice 
spacing in the z direction. Experimentally, we can probe the 
energies of the interacting problem (benzene has 42 electrons) 
by exciting a particle to a reference orbital using a microwave 
transition, which also changes the internal state. This avoids 
problems with postinteraction, assuming that the particles in 
the two different internal atomic states do not interact. Fur¬ 
thermore, the population in this internal atomic state serves 
as an observable measuring the transfer rate as a function of 
excitation energy. 

To demonstrate how interactions change the electronic con¬ 
figuration, let us turn to the conjugated 7 r system, i.e. the 
p z orbitals with delocalized particles, and neglect all interac¬ 
tions with the lower-lying orbitals. When restricting ourselves 
to this subsystem, the problem can be solved using the full 
configuration-interaction method. Experimentally, this sub¬ 
system is also accessible by populating only the six lowest s 
orbitals that are equivalent to the p z orbitals in the x and y 
direction. In the case of p z orbitals, the experimental system 
would incorporate the interaction with the energetically lower- 
lying particles constituting a much more complex situation. 



FIG. 3. Interactions in the conjugated 7 r system, (a) Single-particle 
energy spectrum showing the lowest orbitals (p z orbitals in red) as 
a function of the lattice depth Vo (taking the role of the effective 
nuclei distance Ro). See Appendix for parameters. The energy of 
the p z orbitals can be tuned independently by lattice depth V 0z (here 
35Er). (b) Close-up showing only the p z orbitals with three bond¬ 
ing and three antibonding orbitals, where the orbitals 2-3 and 4-5 
are degenerate, (c) Energy spectrum of the conjugated system incor¬ 
porating the on-site interaction U for Vo = 6Er (see Appendix). 
The colored arrows indicate the lowest excitations (labeled in (e)) 
for vanishing interaction, (d) Logarithmic density plot of the two- 
dimensional Wannier function for the p z (or s ) band. Each contour 
indicates a density drop by a factor 1/10. (e) Level scheme depicting 
orbitals and spin configuration. 


For the description of strongly correlated electron systems, it 
is often desirable to switch to Wannier orbitals (see Appendix) 
that are localized at individual sites (potential minima). The 
Wannier orbital shown in Fig. 3d allows the formulation of 
a tight-binding Hubbard Hamiltonian taking into account the 
tunneling element t\ and the on-site interaction U , where U 
scales linearly with the scattering length a s . The next-nearest- 
neighbor tunneling t 2 is considerably smaller than 1 1 but large 
enough to cause an asymmetry in the p z band (E7 2 ^£o i n 
Fig. 3e). The fermionic many-particle Hamiltonian reads 

H =~E td cl,aC[i+d],<T + c.c. + U^2 (3) 

i,cr,d i 

where the operator c^ a annihilates a particle on site i in spin 
state a G {t,|} and = c\ a Ci j(T corresponds to the re¬ 
spective particle number. Here, the tunneling matrix elements 
t\ (d = 1) and t 2 id = 2) are included and [i] = i mod 6 in¬ 
corporates the periodic boundary conditions. 

For vanishing interaction (a s = 0), the energy spectrum 
(Fig. 3c) represents all possible single-particle excitations of 
the half filled band with delocalized particles (Fig. 3e). In- 
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creasing the interaction causes two drastic changes on the 
’’electronic” structure. First, the particles undergo the tran¬ 
sition from delocalized orbitals to the Mott insulator state, 
where the particles occupy the localized Wannier orbitals. At 
the same time, an interaction gap is formed separating the 
many-particle bands by the on-site interaction energy U. Sec¬ 
ond, the spins of the particles align in an alternating order on 
the Wannier orbitals (antiferromagnetic alignment), thereby 
gaining the second-order energy 2 t\/U per particle. The low¬ 
est Mott band contains possible spin excitations from this 
antiferromagnetic ground state. For the interaction strength 
a s = 0.026, the ratio Ujt\ « 4 is in the vicinity of the pro¬ 
posed values for aromatic hydrocarbons which are still under 
debate [33, 34]. The tunability of interaction strength and po¬ 
tential depth allows an independent control of the ratios V jt\ 
and t^jtx (or the next-neighbor interaction for long-rage in¬ 
teractions [32]). 


IV. CORRELATED ELECTRON DYNAMICS 



One of the outstanding advantages of ultracold atomic sys¬ 
tems is that the dynamical time scale is on the order of milli- 
or microseconds. Therefore, artificial molecules would allow 
monitoring of dynamical processes taking place on femto- or 
attoseconds time scales in real molecules. As an example, 
we compute the hole dynamics after the removal of a particle 
(see Fig. 4). Again, we restrict ourselves to the conjugated 7r 
system, for which the dynamics can be computed. The site- 
selective manipulation has been demonstrated for an optical 
lattice using a tightly focused addressing beam [26, 35]. After 
an evolution time t, suddenly ramping up the optical lattice 
allows freezing the particle number distribution at the indi¬ 
vidual sites of the ring. Using fluorescence imaging, we can 
count the occupation number on each site of the 

artificial molecule [20, 21]. The spin state can be identified 
by selectively removing one of the spin states before imaging 
via resonant light [35]. 

In Fig. 4 the site-resolved dynamics is plotted after the 
removal of a spin-up particle at site vm 3. Without inter¬ 
actions (Fig. 4a), the spin-down particles are not influenced 
by the particle removal, whereas the initial hole in the spin 
up component oscillates between sites i = 3 and 6. This 
quantum revival occurs at the revival time T = H/Eq (with 
Eo = 0.13 Er). This corresponds to the energy difference be¬ 
tween A 2u and Ei g orbitals. For 6 Li, the revival time corre¬ 
sponds to 0.27ms for a recoil energy Er = hx 29.4kHz at a 
lattice wavelength A = 1064 nm. Switching on the (repul¬ 
sive) interaction among the particles (Fig. 4b), the dynamical 
behavior becomes very complex. For t = 0, the antiferromag¬ 
netic order is apparent in the spin-down component causing 
the almost vanishing total density on site i = 3. The neigh¬ 
boring sites are highly occupied with the down component, 
which fills the empty site on short time scales. This inter¬ 


FIG. 4. Correlated dynamics after particle removal, (a) Quantum 
revival for the noninteracting system and (b) dynamics in the antifer- 
romagnetically correlated state for a s /b — 0.02 and Vo = 6Er. The 
plots display the occupation number of the Wannier state at site i of 
the ring for spin up (left) and spin down (right) as a function of time 
t. The images on the left depict the summed density of both spin 
states. The three-dimensional isosurface is plotted for 6/a 2 6 , where 
the coloring refers to the two-dimensional density. 


feres strongly with the dynamics of the up component leading 
to complex density modulations. A priori , it is not clear if 
the spin order persists at the typical temperatures of ultracold 
experiment. Therefore, the dynamics shown in Fig. 4 incor¬ 
porates a finite temperature of 0.05 F/r/&b with k b the Boltz¬ 
mann constant, which corresponds to 70 nK for 6 Li, demon¬ 
strating that the quantum simulation of electronic dynamics is 
feasible with state-of-the-art temperatures. 


V. CONCLUSIONS 

In conclusion, ultracold atoms have great potential as a 
quantum simulator for molecular physics, where solving the 
full many-particle problem is still out of reach using either 
classical or quantum computers [12]. This includes the for¬ 
mation of molecular orbitals, electronic interactions as well 
as the time-resolved monitoring of a broad range of dynam¬ 
ical processes. With the proposed quantum emulation, we 
hope to shed light on open questions regarding the influence of 
electronic correlations and localization in molecular systems. 
In particular, the ability of directly observing spatial corre¬ 
lations goes far beyond the available techniques in molecu- 
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lar physics, where correlations are probed indirectly by pho¬ 
toionization processes. This also touches on the fundamental 
questions of how interactions change the picture of effective 
molecular single-particle orbitals. The electronic correlations 
can fundamentally affect the electron dynamics in a molec¬ 
ular system. While we discuss the creation of an artificial 
benzene molecule using state-of-the-art experimental tech¬ 
niques, the results can be easily generalized covering more 
complex molecules. In particular, spatial light modulators can 
be used to create arbitrary two-dimensional molecular struc¬ 
tures [28]. To emulate nonplanar molecules, one could utilize 
three-dimensional nonseparable lattice potentials formed by 
several interfering laser beams [36] or by holographic projec¬ 
tion techniques. The nuclei positions, kept fixed in the present 
study, can be controlled using slowly varying time-dependent 
optical potentials (e.g. using a spatial light modulator or an ac¬ 
cordion lattice [37]) or coupled to phonons using other atomic 
species such as ions [38]. This allows simulation of the elec¬ 
tronic response to vibrational modes. In principle, one could 
also imagine accessing regimes where the Born-Oppenheimer 
approximation breaks down. 

ACKNOWLEDGMENTS 

We acknowledge funding by the Deutsche Forschungsge- 
meinschaft (Grants SFB 925 and the The Hamburg Centre for 
Ultrafast Imaging (CUI)). 

Appendix A: Proposed experimental sequence 

The optical potential for the artificial benzene molecule is 
generated by three laser beams in the x-y plane, which cre¬ 
ate a honeycomb hexagonal optical lattice Vo Vh C {x, y), and an 
orthogonal one-dimensional lattice Vq z cos 2 ( 7 rz/b) (see Sup¬ 
plemental Material [25]). Superposing this lattice potential 
with a tight dipole trap V& v (x,y) = — CdipVbe -2 ^ +v 
centered on a plaquette of the lattice as sketched in Fig. la 
causes a finite-size system with a hexagonal inner and outer 
ring, where Vo is the height of the potential (in units of the 
recoil energy Er = h 2 /2X 2 m with wavelength A and atomic 
mass m). For concreteness, we choose if not stated other¬ 
wise the lattice depths Vo = 11 Er and Vo z = 35Er, the rel¬ 
ative strength edip = 10 of the dipole potential and the width 
rcdip = 4a of the dipole trap. Here, a — 2A/3 is the width of 
the unit cell in the honeycomb lattice, whereas b = A/2 is the 
lattice spacing in the z direction. The orbital order can be 
tuned by adjusting the lattice depths Vo and Vq z (see Fig. 3a). 

A controlled number of ultracold fermionic atoms can be 
prepared in the tight two-dimensional potential formed by a 
dipole trap and a single antinode of the orthogonal lattice. 
Atoms above a certain energy are expelled from the trap via a 
strong magnetic field gradient [22]. When adiabatically ramp¬ 


ing up the honeycomb lattice, the particles fill up the lowest 
orbitals of the artificial molecule. To image the momentum 
distribution a time-of-flight measurement is performed, where 
the atoms expand for a certain time. Subsequently, a deep 
so-called pinning lattice is ramped up that holds the atoms 
in place while they are illuminated by an optical molasses 
[20, 21]. The pinning lattice can be a square lattice of the 
same wavelength A containing P 2 central lattice sites. Using 
a high-resolution imaging system the individual sites of the 
pinning lattice can be resolved [20, 21]. In addition to the 
observed momenta, the number of particles is also retrieved. 

For observing the momentum density in the z direction, a 
three-dimensional expansion can be performed. In this case, 
the imaging (from the x direction) needs a large depth of field 
and is therefore restricted to a larger optical resolution. How¬ 
ever, for dilute filling of the optical lattice, the atomic posi¬ 
tions can still be identified even if the optical resolution is 
several times the lattice spacing [39] (see Supplemental Ma¬ 
terial [25]). 


Appendix B: Computation of orbitals 

We compute the molecular orbitals by applying the plane- 
wave expansion method for nonperiodic systems. By decom¬ 
posing the two-dimensional potential in Fourier components 
V(x,y) = ^2 v q |q)> the Hamiltonian matrix can be calcu¬ 
lated in the plane-wave basis (2). Here, S is the number of 
unit cells per dimension, q= (. q x ,q y ) and qi = — 1/2+v/S+d 
are the wave vectors with d = {—D, — 79+1,..., D} and 
v = 0,1,..., S-l(i/= 1/2,..., 5-1/2 for S odd). The eigen¬ 
vectors of the Hamiltonian matrix of dimension (2D + 1) 4 5 4 
represent the solution for the single-particle orbitals in the 
plane-wave basis |q). The orthogonal optical lattice in the 
z direction separates. Assuming that only the central lattice 
plane is filled in the experiment, we can restrict ourselves to 
the Wannier function (z) of band n in the z direction. For 
the orbitals in Figs. 1 and 3 we use S = 5, D = 7 and S = 15, 
D — 5 for Fig. 2a. 


Appendix C: Molecular Wannier functions 

The Wannier orbitals at site i of the C ring can be con¬ 
structed from the six s (or p z ) orbitals |V^) via the ba¬ 
sis transformation | Wi) = J2j c ij IV^) with suitable coef¬ 
ficients Cij (see Supplemental Material [25] and Ref. [40]). 
For the lattice depths Vo = 6E R and Vq z = 2bE R (e^p = 10, 
^dip = 4a), the tunneling element is t\ = 0.154^, the next- 
nearest-neighbor tunneling £2 = —0.0094^ and the on-site 
interaction U = 32.9 E R a s /b (U( Pz ) = 22.7 E R a s /b for the p z 
orbital), where a s is the 5-wave scattering length and b = A/2 
(cf. Fig. 3c, d and Fig. 4). 
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